#!/usr/bin/Rscript
cat('Loading packages and data...')
sink('/dev/null')
library(ncdf)
library(methods)
library(spam)
library(fields)
sink()

#get usr requested times/domain
args = commandArgs(trailingOnly=T) #argument = cg#
first.time=as.character(args[1])
last.time=as.character(args[2])

#set grid
north=43; south=39; east=-122; west=-125

#load wind data, remove files from target directory
load('wind.RData') #wind[index]=list(time, lat, lon, u, v)

windDir='../adcirc/'
#system(paste('rm ',windDir,'/*',sep=''))

cat('Done.\n')



#search for requested times
first.time.found = 1
last.time.found = length(wind)
for (t in 1:length(wind)) {
	if (wind[[t]][['time']] < first.time) { first.time.found = t+1 }
	if (wind[[t]][['time']] > last.time) { last.time.found = t-1; break }
}
cat('closest match: begin = ',wind[[first.time.found]][['time']],', end = ',wind[[last.time.found]][['time']],'\n')
if(first.time.found>=last.time.found) { cat('***error in requested times***'); stop() }

#chop wind data
wind=wind[first.time.found:last.time.found]

#find lat/lon of first timestep
lat=wind[[1]][['lat']] #values along first row
lon=wind[[1]][['lon']] #values along first column
u.grid=v.grid=array(dim=c(length(lon),length(lat),length(wind)))
u.grid[,,1]=(wind[[1]][['u']])
v.grid[,,1]=(wind[[1]][['v']])

#generate u/v arrays
for (t in 2:length(wind)) {
	#get lat/lon
	lat.new=wind[[t]][['lat']] #values along first row
	lon.new=wind[[t]][['lon']] #values along first column
	if (any(lat.new!=lat) | any(lon.new!=lon)) {cat('***lat/lon grid changed during time window***'); stop() }
	#get and store u/v arrays
	u.grid[,,t]=(wind[[t]][['u']])
	v.grid[,,t]=(wind[[t]][['v']])
}

#interpolate to adcirc grid
cat('Interpolating to adcirc grid...')
zz <- file(paste(windDir,'/fort.14',sep=''),open='r')
dummy = readLines(zz,n=1)
nnode = scan(zz,nlines=1,quiet=T)[2]
xy = matrix(scan(zz,nlines=nnode,quiet=T),ncol=4,byrow=T)[,2:3]
close(zz)
    
u = matrix(0,nrow=nnode,ncol=dim(u.grid)[3])
v = matrix(0,nrow=nnode,ncol=dim(u.grid)[3])

for(t in 1:length(wind)){
    obj = list(x = lon, y = lat, z = u.grid[,,t])
    u[,t] = interp.surface(obj,xy)
    obj = list(x = lon, y = lat, z = v.grid[,,t])
    v[,t] = interp.surface(obj,xy)
}
cat('Done.\n')

#write the fort.22 file which can be quite large
outfile=paste(windDir,'/fort.22',sep='');
cat('Writitng',outfile,'...\n')

pb <- txtProgressBar(min=1,max=length(wind),style=3)
con <- file(outfile,'w')

for(t in 1:length(wind)){
    
    setTxtProgressBar(pb, t)

    data = cbind(1:nnode,u[,t],v[,t],rep(1013,nnode))  
    data[is.na(data)]=0
    if(t>1)write(paste(' #',wind[[t]][['time']]),file=con,append=T)

	for(j in 1:nrow(data)){
            write(sprintf("%8i%13.5E%13.5E%13.5E",data[j,1],data[j,2],data[j,3],data[j,4]),file=con,ncolumns=4,append=T)
    }

}
close(pb)
close(con)

